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ABSTRACT 

Interstellar dust in galaxies can be traced either through its extinction effects on the star light, 
or through its thermal emission at infrared wavelengths. Recent radiative transfer studies of 
several nearby edge-on galaxies have found an apparent inconsistency in the dust energy bal¬ 
ance; the radiative transfer models that successfully explain the optical extinction underesti¬ 
mate the observed fluxes by an average factor of three. We investigate the dust energy bal¬ 
ance for IC4225 and NGC 5166, two edge-on spiral galaxies observed by the Herschel Space 
Observatory in the frame of the H-ATLAS survey. We start from models which were con¬ 
strained from optical data and extend them to construct the entire spectral energy distribution 
of our galaxies. These predicted values are subsequently compared to the observed far-infrared 
fluxes. We find that including a young stellar population in the modelling is necessary as it 
plays a non-negligible part in the heating of the dust grains. While the modelling approach for 
both galaxies is nearly identical, we find two very different results. As is often seen in other 
edge-on spiral galaxies, the far-infrared emission of our radiative transfer model of IC 4225 
underestimates the observed fluxes by a factor of about three. For NGC 5166 on the other 
hand, we find that both the predicted spectral energy distribution as well as the simulated im¬ 
ages match the observations particularly well. We explore possible reasons for this difference 
and conclude that it is unlikely that one single mechanism is the cause of the dust energy bal¬ 
ance problem in spiral galaxies. We discuss the different approaches that can be considered in 
order to get a conclusive answer on the origin this discrepancy. 

Key words: radiative transfer - dust, extinction - galaxies: structure - galaxies: individual: 
IC4225, NGC 5166 


1 INTRODUCTION 

Although interstellar dust accounts for only a small fraction of the 
interstellar medium (ISM), it plays an important role in various 
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chemical and physical processes in galaxies. Apart from scatter¬ 
ing and absorbing the light at ultraviolet (UV), optical and near- 
infrared (NIR) wavelengths, dust acts as a regulator of the physical 
and chemical conditions in the ISM. It is thus important to under¬ 
stand the amount and distribution of interstellar dust in galaxies. 

When spiral galaxies are seen edge-on, their dust content often 
shows as a distinctive narrow lane where almost all of the stellar 
light is attenuated. Edge-on galaxies allow us to observe both the 
radial and vertical extent of the dust and stellar disk. Because the 
extinction effects are so clearly shown in edge-on spiral galaxies, 
they are the ideal targets when one wants to determine the dust 
content and distribution. 

Studying the dust properties in extinction from UV, optical 
or NIR data offers some clear advantages. The extinction proper¬ 
ties of interstellar dust in this wavelength region are relatively well 
determined and, contrary to the far-infrared (FIR), data with ex¬ 
cellent spatial resolution are easily obtained. The downside of this 
approach is the difficulty to determine the dust mass, as the com¬ 
plex star-dust geometry and physical processes like absorption or 
multiple anisotropic scattering make the conversion from attenua¬ 
tion to dust mass challenging (Disney et al. 1989; Witt et al. 1992; 
Byun et al. 1994; di Bartolomeo et al. 1995; Baes & Dejonghe 
2001; Inoue 2005). In cases where the contrasting light source is 
not the galaxy itself but a known background source, this conver¬ 
sion is more straightforward (White & Keel 1992; Holwerda et al. 
2009, 2013). In the last decade, much progress has been made in 
the field of 3D dust radiative transfer and a number of codes are 
now able to effectively incorporate the necessary physics (e.g., Gor¬ 
don et al. 2001; Baes et al. 2003; Jonsson 2006; Bianchi 2008; Ro- 
bitaille 2011; Schechtman-Rook et al. 2012). An overview of the 
recent developments in 3D dust radiative transfer can be found in 
Steinacker et al. (2013). 

Another way to determine the total dust mass is by measur¬ 
ing the thermal emission of the dust grains at mid-infrared (MIR), 
FIR and submm wavelengths. While there is still some uncertainty 
about the emissivity of dust grains at FIR/submm wavelengths, 
global dust masses can be estimated either using simple modified 
blackbody fits or more advanced SED modelling tools. However, 
in this case, observations are less straightforward and the poor spa¬ 
tial resolution limits detailed information on the dust distribution 
to Local Group galaxies (Meixner et al. 2010; Xilouris et al. 2012; 
Fritz et al. 2012; De Looze et al. 2012b; Smith et al. 2012a) or other 
nearby galaxies such as these from the Very Nearby Galaxy Survey 
(Bendo et al. 2010; Parkin et al. 2012; Bendo et al. 2012; Foyle 
et al. 2012; Hughes et al. 2014). 

As the starlight absorbed by dust grains in the UV/optical/NIR 
region is re-emitted in the MIR/FIR/submm wavelength regime, 
one would expect the absorbed luminosity to be equal to the to¬ 
tal thermal IR/submm luminosity. However, there seems to be a 
deficit in this dust energy balance: the FIR emission predicted from 
extinction in edge-on spiral galaxies seems to underestimate the ob¬ 
served value by a factor of about three (Popescu et al. 2000, 2011; 
Misiriotis et al. 2001; Alton et al. 2004; Dasyra et al. 2005; Bianchi 
2008; Baes et al. 2010; Holwerda et al. 2012b). This is a problem 
typically encountered using self-consistent radiative transfer simu¬ 
lations. SED fitting tools like MAGPHYS (da Cunha et al. 2008) or 
CIGALE (Burgarella et al. 2005; Noll et al. 2009) do not encounter 
this problem. However, these methods cannot take into account the 
effects of the complex star-dust geometry in galaxies, which is nec¬ 
essary for a self-consistent calculation of the dust heating and ther¬ 
mal emission. Radiative transfer codes, on the other hand, can in¬ 
corporate all the geometrical effects, position and strength of the 


heating sources to self-consistently predict the FIR emission. As 
this requires a significant increase in computational cost and with 
the loss of a guaranteed energy balance, they are less useful to study 
larger samples. While methods to predict results of radiative trans¬ 
fer simulations by using artificial neural networks are extremely 
promising, it is advisable to use them for statistical purposes (Silva 
et al. 2012). 

One possible explanation for the discrepancy in the energy 
balance found in radiative transfer simulations is a systematic un¬ 
derestimation of the FIR emission coefficients (Alton et al. 2004; 
Dasyra et al. 2005; MacLachlan et al. 2011). Other authors have 
advocated the presence of additional dust distributed in such a way 
that it hardly contributes to the attenuation, either in the form of 
a second inner disc (Popescu et al. 2000, 2011; Misiriotis et al. 
2001), or as a clumpy or inhomogeneous medium (Bianchi 2008; 
De Looze et al. 2012c, 2014; Saftly et al. 2015). Finally, there is 
now ample evidence for the presence of vertically extended diffuse 
dust in spiral galaxies, which can contribute at least some of the 
infrared emission without contributing much to the large-scale at¬ 
tenuation (Kamphuis et al. 2007; Whaley et al. 2009; Seon et al. 
2014; Hodges-Kluck & Bregman 2014). 

In order to explain this discrepancy, we require a sufficiently 
large set of galaxies that have the needed multi-wavelength data 
available and are modelled in a homogeneous way, such that phys¬ 
ical effects can be disentangled from effects connected to the mod¬ 
elling. Since fitting radiative transfer models to images is a very 
computationally expensive job, most dust energy balance studies 
have concentrated on a single or at most a very modest set of galax¬ 
ies. Moreover, the parameter space of the different galaxy models 
is often degenerate, in the sense that strongly different galaxy mod¬ 
els can result in similar attenuation maps (Bianchi 2007 ; De Geyter 
et al. 2013). Fortunately, these dust energy balance studies can now 
be taken to the next level, thanks to two recent developments. 

The first improvement concerns the determination of the dust 
mass based on optical data. Recently, different attempts have been 
undertaken in the radiative transfer community to make the fitting 
of radiative transfer models to data more automated and objective 
(e.g., Steinacker et al. 2005; Robitaille et al. 2007; Bianchi 2007; 
Schechtman-Rook et al. 2012). One of the challenges here is that, 
in order to fit Monte Carlo radiative transfer models to data, one 
has to take into account that the model images themselves are not 
noise-free, which implies that normal gradient-search optimisation 
methods are problematic. De Geyter et al. (2013) showed that op¬ 
timisation by means of genetic algorithms is an ideal solution for 
this problem. Moreover, the degeneracy in the parameter space can 
be minimised by means of oligochromatic radiative transfer mod¬ 
elling, i.e., by simultaneously fitting radiative transfer models to 
images in several optical wavebands (De Geyter et al. 2014, here¬ 
after DG14). 

The second development is the launch of the Herschel Space 
Observatory (Pilbratt et al. 2010), which has been a major boost for 
the study of the interstellar dust medium in nearby galaxies. Partic¬ 
ularly revolutionary is the SPIRE photometer, which filled a gap in 
the wavelength coverage between previous FIR missions as IRAS 
(Neugebauer et al. 1984), ISO (Kessler et al. 1996) and Spitzer 
(Werner et al. 2004) on the one side, and ground-based submm in¬ 
struments as SCUBA (Holland et al. 1999) and LABOCA (Siringo 
et al. 2009) on the other side. Thanks to Herschel, we can finally 
construct complete dust SEDs, easily trace even the cold dust grains 
and get a more realistic and tighter constraint on the total dust mass. 

The target galaxies in this paper were observed as part of the 
H-ATLAS survey (Eales et al. 2010). This was the largest extra- 
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Property 

unit 

IC4225 

NGC 5166 

“2000 

h 

13:20:01 

13:28:15 

52000 

deg 

31:58:53 

32:01:57 

type 

— 

so/a 

Sb 

Vr 

kms~^ 

5438 

4647 

d 

Mpc 

73.9 

62.6 

scale 

kpc arcsec”^ 

0.358 

0.303 

Dis 

arcmin 

1.3 X 0.3 

2.3 X 0.4 

logZ-K 

^O.K 

10.93 

11.19 

logM* 

Mq 

10.59 

10.82 

log Ltir 

Lq 

9.74 

10.30 

logMd 

Mq 

7.24 

7.89 

log SFR 

Moyr ‘ 

^0.58 

-0.03 

log sSFR 

Gyr-' 

-2.17 

-1.87 


Table 1. The fundamental properties of the two galaxies in our sample. 
Stellar masses, TIR luminosities, dust masses, star formation rates and spe¬ 
cific star formation rates are calculated using a MAGPHYS modelling of 
the spectral energy distribution. 
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galactic survey undertaken by Herschel, with a total of 600 hours 
dedicated observing time. In total, about 600 deg^ of extragalac- 
tic sky was mapped simultaneously with the PACS (Poglitsch et al. 
2010) and SPIRE (Griffin et al. 2010) instruments in 5 bands cen¬ 
tred at 100, 160, 250, 350 and 500 /rm. The scientific goals of the 
H-ATLAS programme are diverse, ranging from the search for de¬ 
bris disks and the characterisation of the Galactic cirrus (Thompson 
et al. 2010; Bracco et al. 2011), statistical studies of the local and 
the high-redshift galaxy population (Dunne et al. 2011; Lapi et al. 
2011; Smith et al. 2012a), to the investigation of extreme objects 
like blazars and gravitational lenses (Negrello et al. 2010; Cox et al. 
2011; Fu et al. 2012; Lopez-Caniego et al. 2013). 

In this paper, we will make use of these two recent devel¬ 
opments to investigate the energy balance in two edge-on spiral 
galaxies, IC 4225 and NGC 5166. In Section 2 we describe the ra¬ 
tionale behind the selection of these two galaxies, and present the 
available data. The construction and results of our radiative transfer 
modelling is presented in Section 3. In Section 4 we discuss these 
results and the possible implications on the dust energy balance 
problem and present our conclusions. 

2 THE TARGETS 
2.1 Target selection 

The starting point for our target selection was the sample of 12 
edge-on spiral galaxies considered by DG14. This sample was se¬ 
lected from the CALIFA survey (Sanchez et al. 2012), based on 
the criteria of an edge-on orientation with a dust-lane morphology, 
a major axis diameter larger than 1 arcmin and a minor axis di¬ 
ameter larger than 8 arcsec. Each of these galaxies were modelled 
in exactly the same way: a stars-l-dust radiative transfer model was 
fitted to the SDSS g, r, i and z images, using the oligochromatic 
FitSKIRT code. 

In this paper, we want to extend this analysis to a panchromatic 
study, and compare the dust “observed” in the optical SDSS images 
to the dust seen in emission in the FIR/submm. We hence select 
from the sample of DG14 those galaxies that have Herschel data to 
cover the FIR/submm SED. 


Figure 1. A comparison of some fundamental properties of the two galax¬ 
ies in our sample and the HRS late-type galaxy population. The top panel 
compares stellar mass and star formation rate, the bottom panel stellar mass 
and dust luminosity. All quantities have been computed in a similar way 
using MAGPHYS. 

Of the 12 galaxies, three galaxies have been observed by Her¬ 
schel: IC2461, IC4225, and NGC 5166. We decided to leave out 
IC2461 from this study - it turned out to be an extreme system 
with a deeply buried Seyfert 2 nucleus (Paris! et al. 2012), that is 
emitting high-energy emission even in the ultra-hard X-ray bands 
(Koss et al. 2011; Baumgartner et al. 2013). Additionally, we lack 
the required UV observations from GALEX, which provides an im¬ 
portant constraint on the younger stellar population. 

Both of the remaining galaxies are located in the North Galac¬ 
tic Plane (NGP) field of the H-ATLAS survey. IC 4225 is situated 
in the Coma supercluster region at a distance of about 73.9 Mpc 
(Slezak et al. 1988; Ramella et al. 1995). It is classified as an SO/a 
galaxy (de Vaucouleurs et al. 1991) and is about 1.3 by 0.3 arcmin. 
NGC 5166 is an Sb type galaxy at a Tully-Fisher based distance 
of 62.6 Mpc (de Vaucouleurs et al. 1991; Theureau et al. 2007). 
It has a prominent but somewhat asymmetric dust lane and has a 
larger angular extent of 2.3 by 0.4 arcmin. Table 1 lists a num¬ 
ber of fundamental properties of the two galaxies. Stellar masses, 
dust luminosities, star formation rates and specific star formation 
rates are calculated by fitting the UV-submm spectral energy dis¬ 
tribution using MAGPHYS (da Cunha et al. 2008). In Figure 1 we 
plot the location of the two galaxies in the M^f^-SFR and the M,),- 
Ltir planes. Overplotted on these panels is the location of the late- 
type galaxy population from the Herschel Reference Survey (HRS, 
Boselli et al. 2010), a complete volume-limited survey of galaxies 
in the local Universe. The plotted quantities are calculated in the 
same way (Viaene et al., in prep.). Table 1 and Figure 1 show that 
IC4225 and NGC 5166 are relatively massive spiral galaxies with 
a high dust content and star-formation rate. This is not unexpected, 
given that one of the selection criteria of the parent sample from 
DG14 was the presence of a clear dust lane. 
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Survey 

wavelength 

[/im] 

IC4225 

[mjy] 

NGC 5166 
[mJy] 

GALEX 

0.226 

0.40 + 0.12 

0.87 + 0.17 

SDSS 

0.354 

1.3+ 0.1 

2.4+ 0.1 

SDSS 

0.475 

5.73 + 0.1 

11.9 + 0.2 

SDSS 

0.622 

12.3 + 0.2 

26.1+0.3 

SDSS 

0.763 

19.0 + 0.3 

40.5 + 0.4 

SDSS 

0.905 

26.6 + 0.3 

55.0 + 0.4 

2MASS 

1.25 

40.3 + 0.5 

80.0 + 0.7 

2MASS 

1.64 

53.0 + 0.5 

127.6+1.0 

2MASS 

2.17 

47.0 + 0.5 

119.6 + 0.9 

WISE 

3.4 

23.0 + 0.2 

51.1 + 0.3 

WISE 

4.6 

13.6 + 0.2 

33.0 + 0.3 

WISE 

12 

15.5 + 0.2 

85.8 + 0.4 

WISE 

22 

11.5 + 0.8 

170 + 3.0 

IRAS 

25 

— 

155 + 38 

IRAS 

60 

— 

845 + 59 

IRAS 

100 

— 

2880 + 259 

PACS 

100 

640 + 72 

3132 + 252 

PACS 

160 

709 + 85 

3506 + 262 

SPIRE 

250 

490 + 35 

2315+ 167 

SPIRE 

350 

190+15 

961+74 

Planck 

350 

— 

1062 + 272 

SPIRE 

500 

54 + 7 

327 + 33 

Planck 

550 

— 

356 + 219 


Table 2. Observed flux densities and corresponding errors for IC 4225 and 
NGC5166. 

2.2 Data 

For the panchromatic radiative transfer modelling that we will per¬ 
form, we need imaging data over a wavelength range as broad as 
possible. Both galaxies have been observed by the Sloan Digital 
Sky Survey (SDSS, Abazajian et al. 2009) and therefore have u, 
g, r, i and z band observations. Additionally they were observed 
by the Galaxy Evolution Explorer (GALEX, Martin et al. 2005), 
the Two Micron All Sky Survey (2MASS, Skrutskie et al. 2006) 
and the Wide-held Infrared Survey Explorer (WISE, Wright et al. 
2010). Unfortunately, neither of the galaxies had a clear detection 
in the GALEX FUV band at 0.152 yum, so only the NUV flux den¬ 
sity of the GALEX survey could be used. For each of these obser¬ 
vations, we followed standard data reduction procedures using IDE 
routines. After masking foreground stars, integrated flux densities 
were determined by means of aperture photometry using DS9 Fun- 
tools (Mandel et al. 2001) instead of using the catalog values to 
ensure that there is no contamination from nearby sources. 

PACS and SPIRE images and flux densities were extracted 
from the H-ATLAS catalogue. Details on the PACS and SPIRE 
data reduction, map-making, source extraction and catalogue gen¬ 
eration for H-ATLAS can be found in Ibar et al. (2010), Pascale 
et al. (2011), Rigby et al. (2011) and Valiante et al. (in prep.). For 
NGC5166 additional data were used from the IRAS Faint Source 
Catalog (Moshir et al. 1992) at 25, 60 and 100 pm. The galaxy was 
also detected at 350 and 550 ^m by Planck where we have used 
the flux density determined using the aperture photometry (APER- 
FLUX) listed in the Planck Catalogue of Compact Sources vl.O 
(Planck Collaboration 2014). 

All the flux densities and their corresponding errors for 
IC 4225 and NGC 5166 are listed Table 2. 


3 PANCHROMATIC RADIATIVE TRANSFER 
MODELLING 

3.1 Modelling approach 

DG14 determined the dust distribution in 12 edge-on spiral galax¬ 
ies, including IC4225 and NGC 5166, based on optical images. 
More specifically, they used the oligochromatic FitSKIRT code to 
construct a radiative transfer model with 19 free parameters that can 
simultaneously reproduce the SDSS g, r, i and z band data for these 
galaxies. For more details on the fitting procedure we refer to Sec¬ 
tion 3 of DG14. The goal of the modelling in the present paper is 
to extend these oligochromatic models to fully panchromatic mod¬ 
els that can reproduce not only the images at optical wavelengths, 
but also explain the thermal emission by dust, as observed in the 
MIR/FIR/submm in a self-consistent way. 

The modelling is done using SKIRT (Baes et al. 2003, 2011; 
Camps & Baes 2015), a Monte Carlo radiative transfer code de¬ 
signed to solve the 3D dust radiative transfer problem. It is mainly 
applied to simulate dusty galaxies (e.g., Gadotti et al. 2010; Gomez 
et al. 2010; De Looze et al. 2012a, 2014), but has also been used to 
model active galactic nuclei (Stalevski et al. 2012) and dusty discs 
around evolved stars (Vidal & Baes 2007; Deschamps et al. 2015). 

Our modelling goes in three phases: our first set of models 
(Case A: Optical-NIR models. Section 3.2) are constrained by the 
optical images only. In a second set of models (Case B: UV-NIR 
models. Section 3.3), we add a component to account for additional 
UV radiation. Finally, we construct a third set of models where we 
add ad hoc components for the cool dust and obscured star forma¬ 
tion in order to fit the excesses that are still present between the 
Case B models and the observed infrared flux densities (Case C: 
UV-NIR ModBB models. Section 3.4). 

3.2 Case A: Optical-NIR models 

The FitSKIRT results from DG14 form the ansatz for the models. 
The stellar distribution consists of an double-exponential disk to 
describe the stellar disk and a flattened Sersic profile to describe the 
bulge. The dust distribution is represented by a double-exponential 
disk with scale length and scale height h, ^. This is a commonly 
used distribution in dust radiative transfer simulations (Xilouris 
et al. 1999; Misiriotis et al. 2001; Tuffs et al. 2004; Bianchi 2007; 
Popescu et al. 2011). Recently, it has been shown, using high res¬ 
olution PACS and SPIRE maps, that such a distribution forms a 
suitable description of the FIR morphology in both radial and ver¬ 
tical direction (Verstappen et al. 2013; Hughes et al. 2014). The 
standard BARE-GR-S model of Zubko et al. (2004) was used to 
set the absorption efficiency, the scattering efficiency and scatter¬ 
ing phase function of the dust. The total dust mass Md completes 
the characterisation of the dust content in the model. 

One remarkable result from this fitting is that IC 4225 has an 
exceptionally large dust disk as can be seen in Table 3. Especially 
when compared in relative size to the stellar disk, it is almost three 
times more extended. It is even more surprising keeping in mind 
that this galaxy was classified as a SO/a galaxy, for which it has been 
shown that the dust emission should be compact (Bendo et al. 2007; 
Munoz-Mateos et al. 2009; Smith et al. 2012b). Most likely how¬ 
ever, due to the smaller angular size, the morphology of IC 4225 
was misclassified. 

The large dust disk, combined with a total dust mass close to 
the average of the modelled sample, results in a face-on V-band op¬ 
tical depth of a mere 0.18. Because of the relatively small angular 
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Figure 2. The SED of IC 4225 (top) and NGC5166 (bottom). The coloured dots with error bars correspond to the flux densities listed in Table 2. The dashed 
lines are the simple models based only on the optical data (Case A), the solid lines are the models in which an additional young stellar component in taken 
into account (Case B), and the dotted lines are the models with additional cold dust and obscured star formation (Case C). The blue transparent band represent 
the spread in the FIR emission of the Case B models within one error bar of the best fitting models. The bottom panels below the SEDs show the residuals 
between the observed SED and the Case B model. 


extent of the galaxy (it is the galaxy with the smallest angular ex¬ 
tent of the sample considered by DG14), the uncertainties on the 
dust parameters are quite large compared to the average values for 
the sample. In spite of the large uncertainties, the model reproduces 
the optical images to an acceptable level with only a small irregu¬ 
larity in the outer regions, possibly due to the projection of a spiral 
structure (see Figure 3; of DG14). 

The model for NGC5166 is more in line with those found 
for other edge-on spiral galaxies. The dust disk is about 50% more 


extended and about half as thick as the stellar disk, which corre¬ 
sponds to the classical image of edge-on spiral galaxies (Xilouris 
et al. 1999; Bianchi 2007). The dust parameters were constrained 
with a higher accuracy compared to IC 4225. The model has a bulge 
that is quite large, relative to the disk. As was found by Liitticke 
et al. (2000b) and can be clearly seen from the residuals in Figure 
3j in DG14, the bulge has a strong peanut shape, a feature that is 
quite common when looking at edge-on spiral galaxies (Liitticke 
et al. 2000a). As was the case for IC4225, a small asymmetry can 
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Parameter 

unit 

IC4225 

NGC 5166 

hR,ii:,ev 

kpc 

3.38 + 0.05 

3.92 + 0.17 


kpc 

0.84 + 0.04 

0.66 + 0.04 

ReS 

kpc 

1.25 + 0.04 

3.17 + 0.59 

n 

— 

1.1 +0.1 

5.8+ 0.9 

4 

— 

0.69 + 0.03 

0.26 + 0.04 

hR4 

kpc 

10.01 + 2.74 

5.81 + 1.12 

hz.i 

kpc 

0.24 + 0.02 

0.34 + 0.06 

Mi 

10^ Mq 

2.1 +0.5 

4.8+ 0.9 

■'"v 

— 

0.18 + 0.17 

1.23+0.37 


— 

7.4+ 0.8 

21.2 + 4.7 

i 

deg 

89.3 + 0.1 

87.6 + 0.5 

1 

IO’Lq 

4.16 + 0.1 

11.39 + 0.51 

BITi 

— 

0.43 + 0.01 

0.45 + 0.02 

SSPdisk 

Gyr 

7 

4 

SSPbulge 

Gyr 

12 

11 

young 

kpc 

10.01 

5.81 

young 

kpc 

0.12 

0.17 

Td.ModBB 

K 

21 

20 

Md.ModBB 

10’ Mq 

0.9 

1.5 


Table 3. The parameters of the best fitting models for 1C 4225 and 
NGC5166. The top 13 parameters with error bars correspond to the Fit- 
SKIRT models as described in DG14. The bottom six parameters are the 
additional parameters introduced for the panchromatic radiative transfer 
models: the SSP ages for disk and bulge for the Case A models, the scale 
length and scale height of the young stellar disc for the Case B models, and 
the parameters of the modified black body dust emission for the Case C 
models. 

be seen, possibly due to a spiral structure. In general, the fit is very 
satisfactory with most pixels only showing a maximum deviation 
of 25%. 

The first step in our analysis consists of adjusting these 
oligochromatic FitSKIRT models in order to calculate the view of 
the galaxy in the entire UV to submm domain, rather than only in a 
few selected optical bands. This implies that both the properties for 
the stars and dust need to be set over this entire wavelength domain. 

For both stellar components, we assume an intrinsic Bruzual 
& Chariot (2003) single stellar population SED, which is based on 
Chabrier (2003) initial mass function and assuming a solar metal- 
licity. The ages of the best fitting SSPs were determined from the 
dereddened luminosities reported in DG14 using an adapted ver¬ 
sion of the SED fitting tool from Hatziminaoglou et al. (2008, 
2009). For IC 4225, we found ages of 7 Gyr for the stellar disk and 
12 Gyr for the bulge component. The best fitting stellar populations 
for NGC 5166 are slightly younger at 4 Gyr for the stellar disk and 
11 Gyr for the bulge component. The entire SED was normalised 
to match the i-band luminosity. 

Concerning the dust, the standard BARE-GR-S dust model 
was used to set the optical properties in the entire UV-submm 
wavelength domain. Note that the SKIRT code considers individ¬ 
ual silicate, graphite and PAH grain populations for the calculation 
of the dust emission spectrum, and that it fully takes into account 
the transient heating of PAHs and very small grains, using tech¬ 
niques described in Guhathakurta & Draine (1989) and Draine & 
Li (2001). 

In Eigure 2 we compare the SED obtained for the Case A 
model (dashed lines) with the observed flux densities for both 
galaxies. It is clear that for both galaxies the models cannot ex¬ 
plain the observed SEDs. The radiative transfer models reproduce 


the optical SED very well, which is not surprising given that the 
geometrical parameters of the model and the characteristics of the 
intrinsic SEDs were determined to reproduce the optical images. 
They can also successfully explain the NIR data points. The model 
SEDs, however, underestimate the observed flux densities both in 
the UV and in the thermal MIR/EIR/submm regions. 

3.3 Case B: UV-NIR models 

It is clear that the simple Case A models considered above need to 
be extended with an additional component that can account for the 
observed UV radiation. We will account for this by adding a third 
stellar component associated with young, unobscured star forming 
regions. Unfortunately, for both galaxies, the GALEX NUV images 
are too noisy to be used to constrain the geometry of the younger 
component without any prior assumptions. Young stellar popula¬ 
tions in spiral galaxies are usually located in a thin disk with scale 
heights similar or smaller than the dust disk (Schechtman-Rook & 
Bershady 2013, 2014). In our models, the scale height of this young 
stellar disk is taken to be half that of the dust disk'. The values we 
use are 120 and 170 pc for IC4225 and NGC5166 respectively, 
similar to what was found as the scale height of the young stellar 
disk in the Milky Way (Bahcall & Soneira 1980). 

For the intrinsic SED of this young component, we used a 
Starburst99 SED template which represents a stellar population 
with a constant, continuous star formation rate and an evolution 
up to 100 Myr (Leitherer et al. 1999). In this case, the initial mass 
function is a Salpeter (1955) IME with masses between 1 Mq and 
100 Mq and again assuming a solar metallicity. The luminosity of 
this component is constrained by the GALEX NUV flux density. 
As the contribution of this SED is negligible in the optical bands, 
adding this component does not affect the models in the SDSS and 
2MASS bands. As it does, however, affect the energy balance of 
the dust, a new set of radiative transfer simulations has to be mn in 
order to constrain the NUV flux and determine the resulting images 
over the entire wavelength range. 

The solid line in the left panel of Eigure 2 compares the model 
SED of this new Case B model with the observed flux densities for 
IC 4225. Compared to the simple model without young stellar com¬ 
ponent, this model is a clear improvement. The optical/NIR SED is 
still reproduced as accurately as in the simple model, but the model 
now also reproduces the UV portion of the SED. As a result of the 
additional energy input of the young stellar populations, much of 
which is absorbed by the dust, the infrared part of the SED is also 
boosted. However, it is clear that a EIR excess still remains: the 
model EIR/submm flux densities still underestimate the observed 
data points by a factor of a few. 

Interestingly, we obtain a different picture for NGC 5166. 
There was also a clear EIR/submm inconsistency between the data 
and the simple Case A model. Contrary to IC4225, however, this 
inconsistency disappears almost completely when we include a 
young stellar component: not only does the model now reproduce 
the UV part of the SED, but it almost provides sufficiently addi¬ 
tional dust emission to account for the observed EIR/submm emis¬ 
sion. The only data points that are still strongly underestimated by 
the model are the WISE 22 /jm and IRAS 25 yum, which are usually 

' For both galaxies we ran additional simulations with different values for 
the scale height of the young stellar' population, ranging between one third 
of the dust scale height to the same value. This effect, however, was hardly 
noticeable, as was also found by De Looze et al. (2014). 
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linked to obscured star formation (Calzetti et al. 2007, 2010; Lee 
etal. 2013). 

To investigate the effect of the uncertainty on the dust distribu¬ 
tion of IC4225, we have modelled the spread in the FIR emission 
when starting from a model with slightly altered parameter values. 
Table 3 in DG14 mentions both the best fitting values and the un¬ 
certainties for every free parameter. The highest FIR emission can 
be expected to be the model with the smallest dust disk and highest 
dust mass while the lowest FIR emission is found in the opposite 
case. Figure 2 shows the spread on the Case B model where the lu¬ 
minosity of the young stellar component was adjusted to reproduce 
the original value. By using this approach we investigate the FIR 
emission of models that still appropriately model the optical bands 
and make no additional assumptions using MIR/FIR data. 

It is clear that for IC 4225, the spread on the FIR emission is 
larger, resulting from the higher uncertainty on the dust distribu¬ 
tion. The observed data points are only barely compatible with the 
densest dust disk model. However, in this case even the 22 fim flux 
is recovered. For NGC 5166, the spread in the FIR is considerably 
smaller and not much variation on the shape of the SED can be 
expected. 

In Figures 3 and 4 we compare a set of images correspond¬ 
ing to the Case B models to the actual observed images of the two 
galaxies in different wavebands across the entire wavelength range 
calculated in the simulations. The simulated images have been con¬ 
volved with the appropriate point spread function: for GALEX, 
SDSS, 2MASS and WISE we used a simple Gaussian with the cor¬ 
rect FMHM, for the PACS and SPIRE data we used the kernels 
provided by Aniano et al. (2011). 

From the GALEX NUV image it is clear that for IC 4225, the 
unobscured young stellar component is considerably more concen¬ 
trated in the center, than what is assumed in our model. Addition¬ 
ally, the MIR and FIR images seem to suggest that, especially the 
warm dust, is also concentrated closer to the centre of the galaxy. 
It is clear that the dust disk as obtained in DG14, is clearly not a 
realistic model for this galaxy. While it successfully reproduces the 
optical images, the disk is too extended to reproduce the UV and 
MIR/FIR observations. Interestingly, our radiative transfer model 
was capable of recovering the 500 ^m flux density. However, it 
seems that the model does not really predict the correct morphol¬ 
ogy: the model image at 500 /rm is spatially extended, while the 
observed SPIRE 500 /rm image seems much more centrally con¬ 
centrated. 

For NGC 5166 the comparison between model and observa¬ 
tions is much more satisfactory. In general, the geometry of the 
model images seems to match the observations reasonably well. 
While less obvious due to the low signal-to-noise image, the NUV 
images seem to be in better agreement than was the case for 
IC 4225. Apart from the boxy/peanut shaped bulge and some small 
bright, additional features, probably projected spiral arms, the near- 
infrared images also match quite well. These off-centre features 
also seem to be slightly noticeable at 100 ^m but, in general, the 
model is in good agreement with the observations. 

3.4 Case C: UV-NIR + ModBB models 

To get some more insight into the missing component in our model 
we have added two additional components to the previous model. 

First, to explain the underestimated emission in the 
FIR/submm regime, we fit a modified black body through the resid¬ 
ual Herschel fluxes and add this component to our model. As in 
both galaxies the 500 /rra flux was recovered, this leads to a small 


overestimation of the cold dust contribution in this model. From 
this fit we can determine the amount of cool dust that is missing 
and the associated dust temperature. For IC 4225 we require an ad¬ 
ditional 9x10*’ Mq at a temperature of 21 K to explain all FIR data 
points, which amounts to about one third of the total dust mass. A 
significant part of the dust of IC 4225 should hence be distributed in 
such a way that it does not add significantly to the large-scale atten¬ 
uation but has a significant contribution in the FIR part of the SED, 
such as quiescent compact clumps. Combining both components, 
the total dust mass of NGC 5166 is 3.0 x 10’ Mq, 75% larger than 
the value estimated using MAGPHYS (1.7 x lO’ Mq, see Table 1). 
The difference is probably the result of the overabundance of cold 
dust in our model (due to the extended dust disk, our model over¬ 
estimates the observed flux densities at the longest wavelengths). 

In the case of NGC 5166, the best fit returns 1.5 x lO’ Mq of 
dust at a temperate of 20 K. However, it should be noted that this 
overestimates some of the observed Herschel fluxes and should be 
considered as an upper boundary for the missing dust mass. There¬ 
fore, at most one fourth of the total dust mass of NGC 5166 is 
locked up in quiescent compact clumps. In this case, the total dust 
mass amounts to 6.3 x 10’ Mq, in very good agreement with the 
MAGPHYS value of 7.7 x 10’ Mq (Table 1). 

Secondly, the WISE 22 /rm and IRAS 25 yum mid-infrared 
fluxes of both galaxies are underestimated by the model. Emis¬ 
sion at these wavelengths is dominated by obscured star formation. 
Some obscured star formation is already included in the Case B 
models. Indeed, the observed NUV flux was fitted by adding a 
young stellar component, and a significant part of the correspond¬ 
ing UV emission was attenuated by the diffuse dust component. 
This additional dust heating and corresponding mid-infrared emis¬ 
sion is clearly visible as the difference between the Case A and 
Case B lines in Eigure 2, but not sufficient to explain the ob¬ 
served mid-infrared flux densities. The remaining mid-infrared ex¬ 
cess can be modelled as completely obscured star forming regions 
(De Looze et al. 2014). We could opt to add another modified black 
body, but we preferred to use an ultra-compact Hn-region template 
from Groves et al. (2006). Since this component is only constrained 
by the 22 and 25 yum fluxes, the corresponding dust masses are 
poorly constrained (in any case, they are negligible compared to 
the cool dust mass budget). 

The Case C models for both galaxies are represented as the 
dotted lines in Figure 2. We would like to stress that this model is 
an ad hoc model that holds no “predictive” power. This is unlike 
previous models, where the entire SED is essentially derived in a 
self-consistent way from NUV/optical information only. 


4 DISCUSSION AND CONCLUSIONS 

The ambition of this paper was to use the state-of-the-art, both ob- 
servationally and in radiative transfer modelling, to investigate the 
so-called dust energy balance problem for edge-on spiral galax¬ 
ies. Eor both galaxies, IC4225 and NGC 5166, the observed SED 
shows a clear FIR excess compared to a standard model that is 
based only on optical images, and hence only contains evolved 
stellar populations. This is similar as what is found by various 
teams (e.g., Popescu et al. 2000; Misiriotis et al. 2001; Alton et al. 
2004; Dasyra et al. 2005; Bianchi 2008; Baes et al. 2010; Popescu 
et al. 2011; Holwerda et al. 2012b). However, things change when 
a young stellar population is included in the modelling: while a 
substantial FIR excess still remains for IC 4225, the FIR fluxes are 
recovered almost exactly for our model of NGC 5166. For the lat- 
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Figure 3. Comparison between the observations (left) and model images (right) for IC4225. All images correspond to the Case B model 
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Figure 4. Same as Figure 3, but now for NGC 5166. 


ter galaxy, the Case B model does not only reproduce the observed 
SED, the predicted images over the entire UV-submm wavelength 
range also provide a good match to the observed images. For both 
galaxies, the SED can be forced to fit the data by adding an ob¬ 
scured star formation template and a component of cool dust that 
dust not affect the large-scale extinction. The latter could be inter¬ 
preted as a population of dense and quiescent clumps. 

Why these two galaxies show such a different behaviour is an 
interesting question. It is possible that this is a consequence of the 
modelling. For both galaxies, we used exactly the same approach, 
adopted the same assumptions and used the same data as input. 
Moreover, the FitSKIRT fitting of radiative transfer models to the 
optical images in both cases resulted in fits of similar (and satis¬ 
factory) quality. One aspect that is different, though, is the different 
angular size of the galaxies. 

IC4225 is the smallest galaxy, in angular size, of the sam¬ 
ple studied in DG14. As a result, some of the parameters in the 
model were not constrained as well as for most of the other galax¬ 
ies. In particular, the derived intrinsic properties of the dust disk in 
IC 4225 are rather exceptional: the dust disk scale length retrieved 
is approximately 10 kpc, which is especially large when compared 
to the more moderate stellar disk with a scale length of 3.4 kpc. 
However, when using the uncertainties on the derived parameters 
from optical modelling, the model with the highest FIR output is 
just within range of the observed values. 


NGC 5166, on the other hand, is one of the largest galaxies in 
the sample, and the recovered star-dust geometry is more in line 
with the average results of other studies of edge-on spiral galaxies 
(Xilouris et al. 1999; Bianchi 2007). While this difference in ap¬ 
parent size and hence resolution might be partly responsible for the 
different results, it remains strange that NGC 5166 shows hardly 
any FIR/submm excess. Indeed, this is the galaxy with the most 
“regular” star-dust geometry, and one would expect that this galaxy 
would show a FIR/submm excess as most other “regular” edge-on 
spiral galaxies. The main difference between our analysis, com¬ 
pared to previous work on other edge-on galaxies, is a more self- 
consistent radiative transfer modelling in the optical, since we fit 
our models to multiple bands simultaneously. Whether or not this 
is an important aspect remains to be investigated. 

While we have modelled both galaxies in the same way, it 
is not impossible that these two systems are intrinsically different. 
For example, the compactness of the UV and FIR emission in both 
galaxies is fairly different: for IC 4225 the emission is considerably 
more centrally concentrated than for NGC 5166. Also the spatial 
distribution of the stars is not the same: the Sersic indices of the 
bulges are noticeably different, at 1.1±0.1 and 5.8±0.9 for IC 4225 
and NGC 5166 respectively. This could point to different evolution¬ 
ary mechanisms of both galaxies that could result in a different ISM 
structure (i.e., with a different degree of dumpiness and substruc- 
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ture), and a difference in the ratio of obscured versus unobscured 
star formation. 

The main conclusion of our investigation is that, instead of 
gaining more insight in the driving mechanism for the dust energy 
balance problem, we seem to find that the problem might be more 
complex than anticipated or hoped for. If we want to get a more 
detailed insight in the mechanisms involved in the dust energy bal¬ 
ance, it is necessary to perform this panchromatic radiative transfer 
modelling on a larger scale. Up to now, such modelling has almost 
exclusively been applied to individual galaxies, with different mod¬ 
els, assumptions and codes making the comparison or search for 
correlations essentially impossible. A uniform modelling effort on 
a larger sample of edge-on galaxies would be highly desirable. 

One necessary tool for such an effort is the availability of ra¬ 
diative transfer fitting codes that can efficiently search large pa¬ 
rameter spaces and work semi-automatically; the work by DG14 
demonstrates that this is nowadays possible. Another crucial ingre¬ 
dient is obviously the availability of suitable targets. Excellent cat¬ 
alogs of edge-on spiral galaxies are available nowadays (Kautsch 
et al. 2006; Bizyaev et al. 2014), and the fact that relatively large 
areas of the extragalactic sky have been observed with Herschel in 
the frame of projects such as H-ATLAS opens possibilities. 

The combined HEROES (Verstappen et al. 2013) and 
NHEMESES (Holwerda et al. 2012a) samples contain 19 nearby 
edge-on spiral galaxies that have been observed with PACS and 
SPIRE. These galaxies are the ideal candidates for a follow-up 
study as they have a spatial resolution that is an order of magnitude 
better than the two galaxies studied in this paper. We are currently 
collecting the necessary optical data to start the modelling of this 
sample. 

Apart from systematically modelling larger sets of edge-on 
spiral galaxies, there are other options that need further exploration. 
One interesting option is looking at face-on spiral galaxies (De 
Looze et al. 2014). In this case, smaller scale structures remain 
present in the observations instead of getting smoothed out over 
the line of sight. The disadvantage of this approach is that there is 
no direct analogue for the dust lanes in edge-on spiral galaxies, by 
which the dust mass can be estimated. Additionally, because of the 
lack of information on the vertical distribution, assumptions have to 
be made. These can be extracted from systematic studies of edge-on 
spiral galaxies (Xilouris et al. 1999; Bianchi 2007; De Geyter et al. 
2014; Bizyaev et al. 2014), with additional priors or constraints 
provided by occulting galaxy studies (White & Keel 1992; Holw¬ 
erda et al. 2009, 2013; Holwerda & Keel 2013) or power spectmm 
analyses (Block et al. 2010; Combes et al. 2012). 

Another necessary follow up study is the influence of the as¬ 
sumed dust model. While it is clear that changing the dust opaci¬ 
ties and grain composition might not change the difference between 
both galaxies extremely, the overall influence of the made assump¬ 
tions on the global SED was not investigated in this paper. 

Einally, a different path that needs further inspection is mod¬ 
elling of mock observations using radiative transfer simulated 
galaxies. The use of artificial data, in which the intrinsic proper¬ 
ties are known and compared to those obtained after the radiative 
transfer modelling, can yield interesting clues on biases or system¬ 
atic effects in the modelling. The pioneering work by Misiriotis 
et al. (2000) and Misiriotis & Bianchi (2002), who investigated the 
effects spiral arms and a clumpy ISM respectively, was still based 
on relatively uncomplicated toy models. 

A more advanced option is using the results of 3D hydrody- 
namical simulations, which have managed to produce more real¬ 
istic spiral galaxies in the past few years (Guedes et al. 2011; Re- 


naud et al. 2013; Marinacci et al. 2014; Schaye et al. 2015). Several 
teams have started using such hydrodynamical simulations to inves¬ 
tigate the systematic effects related to radiative transfer effects (e.g., 
Hayward & Smith 2015). Our initial results based on applying Fit- 
SKIRT to two of these models seems to suggest that the complex 
and inhomogeneous dusty medium in spiral galaxies could “hide” 
significant amounts of dust (Saftly et al. 2015). One aspect, very 
relevant for the present study, that could be also investigated using 
such models is the role of resolution in radiative transfer modelling. 

Given the various strengths and weaknesses of each these ap¬ 
proaches, pursuing all three options might be the best approach to 
gain a conclusive answer on the driving mechanisms causing the 
dust energy balance. 
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